Index theoretic characterization of d-wave superconductors in the vortex state 
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We employ index theoretic methods to study analytically the low energy spectrum of a lattice 
d-wave superconductor in the vortex lattice state. This allows us to compare singly quantized hc/2e 
and doubly quantized hc/e vortices, the first of which must always be accompanied by Z2 branch 
cuts. For an inversion symmetric vortex lattice and in the presence of particle-hole symmetry we 
prove an index theorem that imposes a lower bound on the number of zero energy modes. Generic 
cases are constructed in which this bound exceeds the number of zero modes of an equivalent 
lattice of doubly quantized vortices, despite the identical point group symmetries. The quasiparticle 
spectrum around the zero modes is doubly degenerate and exhibits a Dirac-like dispersion, with 
velocities that become universal functions of Ao/t in the limit of low magnetic field. For weak 
particle-hole symmetry breaking, the gapped state can be characterized by a topological quantum 
number, related to spin Hall conductivity, which generally differs in the cases of the hc/2e and hc/e 
vortex lattices. 



Understanding the interactions among fermionic quasi- 
particles and topological excitations of two dimensional 
superconductors constitutes one of the main challenges in 
modern condensed matter physics. It is almost certainly 
true that at least some of the peculiar phenomena ob- 
served in the underdoped cuprates are related to the ex- 
istence of nodal fermionic excitations of a d x 2_ y 2wa.ve su- 
perconductor with low superfluid density, which is even- 
tually susceptible to dephasing by proliferation of mobile 
vortices pj. 

When an external magnetic field induces a finite den- 
sity of vortices which form an Abrikosov lattice, the mo- 
tion of vortices can be neglected and the problem of the 
quasiparticle spectrum is simplified. The initial theoret- 
ical analysis was based on numerical computations Q, 
semiclassical approximation |2|, and scaling arguments 
0. As pointed out by Franz and Tesanovicj^, E3> > the 
quasiparticles interact with the hc/2e vortices by both 
the (semiclassical) coupling of the charge currents to the 
superflow as well as a purely quantum coupling of the 
spin currents to the Zi branch cuts emanating from the 
vortices. Marinelli, Halperin, and Simon [8j then showed 
that within the linearized approximation the spectrum at 
each node remains gapless to all orders of perturbation 
theory in phase gradients. This result was rederived by 
an insightful symmetry analysis in Ref. ||, where it was 
concluded that the effect of the hc/2e vortices is merely a 
possible renormalization of the nodal quasiparticle veloc- 
ities, but that the presence and particularly the number 
of the zero modes is the same as in the lattice of doubly 
quantized hc/e vortices. Various topological aspects of 
the problem were discussed in |7|, |3, [l0( . 

In this Letter, we revisit the problem of the Fermionic 
(quasiparticle) spectrum of a lattice d-wave superconduc- 
tor in the vortex state. Our main result is that there 
can arise a fundamental topological difference between 
the singly and doubly quantized vortex lattice which is 



beyond perturbative analysis [a l2| • ^ n the case of an in- 
version symmetric vortex lattice, this difference relates 
to the non-trivial transformation properties of the Z2 
branchcuts under inversion, absent for hc/e vortices. 

Our point of departure is the inversion symmetric vor- 
tex lattice state of the <i-wave superconductor on a tight- 
binding lattice with perfect particle-hole (p-h) symmetry. 
We use the p-h symmetry to conduct non-perturbative 
analysis of the low energy spectrum and find conditions 
sufficient for the existence of zero energy nodal states. 
Armed with this analysis, we study weak breaking of the 
p-h symmetry, which gaps the nodal points, and show 
that the resulting state is topologically non-trivial. It is 
characterized by a topological quantum number propor- 
tional to the spin-Hall conductivity 0, EJ EH E3 • 

Remarkably, in cases when the branch cuts transform 
non-trivially under inversion, the number of the zero 
modes is doubled for hc/2e vortices relative to the hc/e 
counterpart, despite identical discrete translational sym- 
metry of both vortex lattices in question. The difference 
also holds for a lattice of hc/2e vortices, when the branch 
cuts are ignored. We demonstrate this by proving a form 
of an index theorem [T3L Il4| , which puts a lower bound 
on the number of zero modes, and by explicit analyti- 
cal and numerical solutions, which reveal that the lower 
bound is, in fact, saturated. 

We now provide justification for the above claims. 
Consider the Hamiltonian H = Hq — [iN, where 



Hn = 



(TV') 
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where the sum is over the nearest neighbors (rr') of the 
(underlying) tight-binding lattice. The hopping integral 
i rr ' = t exp(— iA rr r); in the symmetric gauge the mag- 
netic flux $ enters the Peierls factor via A rr +£ = y$/2 
and A rr +y = — x$/2, and the d-wave pairing field in 



the vortex lattice state reads A rr / = ?7 r _ r /Ao exp(i# rr /); 
T]8 = +(— ) if S || x(y). The bond phase factors exp (i6 Tr >) 
result from a self-consistent calculation, but the main 
topological feature of 9 rr i is its 2tt winding around the 
magnetic field induced vortices. As the initial Ansatz for 
the self-consistent solution, we choose the lattice bond 
variables 8 Tr > as follows: first we solve the continuum 
equations VxV^(r) = 2nz ^\ S(r— r^) and V-V0(r) = 
where denotes the vortex positions which in this Letter 
form Abrikosov lattice. These conditions determine </>(r) 
up to (fio + vo ■ r, which is fixed by requiring zero over- 
all current. Next, we define the site variables e l ^ r as the 
value of exp(i</>(r)) at each site, and the initial bond vari- 
able e 10 "' is taken to be the geometric mean of two neigh- 
boring site variables: e ie **' = (e 1 ^ + e l ^')j\e 1 ^ + e^»' |. 
This choice guarantees that if the vortices in e l ^ r reside 
inside plaquettes, so will the vortices in e l9rT ' . 

The diagonalization of the Hamiltonian Q is equiva- 
lent to the solution of the Bogoliubov-de Gennes (BdG) 
equation Ho^pr = Etp r where the lattice operator 



Ho = 



£ r - n 
At 
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Both £ r and A r are defined through their ac- 
tion on a wavefunction at the lattice site r as 



= -*£,= 
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s Vs u r+s- The BdG Hamiltonian Ho, as well 



A £, 

as Pauli a matrices used later, act on the two component 
Nambu spinor ip r = [u r , v r ] T . 

Translational symmetry: While the vortex positions 
are periodic, the Hamiltonian (J2J is invariant only if the 
discrete translations are followed by a gauge transforma- 
tion (magnetic translations). However, as shown in Rcfs. 
B Hi Ho can be transformed into a periodic Hamiltonian 
by a singular gauge transformation, U = cxp{|cx3cV}. 
On a lattice, the phase factors e ± ^ r must be defined 
with care. First, connect pairs of the singly quantized 
vortices by branch cuts (see Fig.l). Then, choose a ref- 
erence point ro at which bo is either one of the two so- 
lutions to 6q = e'^'o, and set e^ r o = b - Now, let r 
be a site neighboring ro not connected by a bond that 
crosses a branch cut, and determine e^ r by the solu- 
tion to b 2 = e l ^ r closer to e^o, i.e. minimizing |6o — b\. 
Next, set r as a reference point and determine the phase 
of its neighbors in the same way. This procedure defines 
e 3^ r on all sites of the lattice uniquely. It also ensures 
that the difference of the phases of e^ r and e^ r ' lies 
inside (—§,§) mod 2tt on all bonds, except across the 
bonds crossing the branch cut where it is more than 
It follows that 



\pi4>r -L pi<t> r 



~2.: 



where Z2,rr' = 1 on each bond except the ones cross- 
ing the branch cut where Z2. rr' = 

— 1. The transformed 
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FIG. 1: Magnetic unit cell (for £—6) containing two hc/2e vor- 
tices represented by two solid (red) circles. The tight-binding 
lattice is shown by the thin (black) lines, the (Z2) branch-cut 
connecting the vortices is shown by a dashed (blue) line. The 
inversion of the branch cut about the midpoint between the 
vortices defines the two regions C (black dots) and C* (white 
dots). The gauge factor 7(1") defined in the text, is 1 for the 
points in C and —1 for the points in C* . 



Hamiltonian, Ho =IA 1 Ho U, is now 

Ho = 0-3 (£r - m) + o"i A r , (3) 
where the transformed lattice operators satisfy 

Srlpr = ~t Yl Z 2,rr+5 X e^ y "-+^ r+5 (4) 
<5=±x,±y 

Ari = A ^ Z 2 ,rr+S X T] S ^ r+S . (5) 
5=±x,±y 

The physical superfluid velocity enters via the factor 



l+e' 1 



|l+e 



—77 



-e iArr ' and describes the lattice ver- 
sion of the semiclassical (Doppler) effect. On the other 
hand, the effect of the Z2 field is purely quantum mechan- 
ical and tied to the hc/2e flux quantization. It would be 
absent for doubly quantized hc/e vortices. 

The resulting Hamiltonian Ho is invariant under dis- 
crete translations by £ x ^ y defining the magnetic unit cell, 
due to the periodicity of V Tr > 0| and our periodic choice 
of the branch cuts. Thus, it can be diagonalized in the 
Bloch basis. Extracting the crystal wavevector k from 
the Bloch wavefunctions gives H(k) = e~ jkr 7Yoe jkr act- 
ing on the Hilbert space of periodic Nambu spinors. 

Inversion symmetry: Our subsequent analysis will fo- 
cus on H(k). For a unit cell with l x l y sites, there are 
2£ x £ y eigenvalues for each k in the first Brillouin zone. 
Consider now an inversion symmetric vortex lattice for 
which e l ^ r = e l ^- r . The inversion operator I, defined by 
its action on the wavefunctions 2/0 r = ip- r , transforms 
the Hamiltonian according to I7i(k)X = H'(— k). The 
important point is that H'{— k) does not always equal 
to H{— k), because the branch-cut can transform in a 
non-trivial way under inversion. A square vortex lat- 
tice depicted in Fig. 1 represents such a case. How- 
ever, the branch cut can be restored by an additional 
Zi gauge transformation j r , where j r — 1 Vr G C, and 



hc/2e G* r = (0,0) |(l, 0) 1(1,1) 1(0, I) 
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(0,0) 
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-1 
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X: 


(-27T/A0) 
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M: 


{-2n/t, -2tt/£) 
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-1 


-1 


-1 
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(0, -&r/*) 
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-1 



hc/e G* r = (0,0) 1(1,0) 1(1,1) 1(0,1) 



r : 


(0,0) 
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-1 
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-1 


X: 


(-2tt/£,0) 
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-1 


-1 


M: 


{-2ir/£, -2-k/I) 


1 


1 


1 


1 


Y: 


(0, -2tt/^) 


1 


-1 


-1 


1 



TABLE I: The value of the diagonal elements of Vu* — 
(— l) x+y e lG r 7(r)Tl at invariant points for the case in Fig^ 
For hc/2e vortices, 7 r is non-trivial and TrPk* = ±4 for all k* 
indicated in the leftmost column. This guarantees at least 16 
zeros in the Brillouin zone. On the other hand, for hc/e vor- 
tices 7 r = 1 and TtV = 8 only at k* = M, which guarantees 
only 8 zeros. For band notation see Fig|5| 

7 r = — 1 Vr e C*. The set C is defined by the points in- 
side the area enclosed by the branch cut and its inverted 
image, and C* is its complement (see Fig. QJ. Thus, 
7 r X7i(k)X7 r = H(— k). Moreover, the spectrum of W(k) 
is symmetric under E —> —E for any /j, |7J. 
Index theorem: For fi = we can define 

T> = (-l)*+" 7p Z 1, (6) 

where 1 is the unity in the Nambu space and the function 
(— l) x+ y changes sign on every other site. Clearly, V 2 
1 . In addition this operator satisfies 

VH(k)V = -H(-k), (7) 

and, in particular, at k = 0, VH(0)V = -H(Q). There- 
fore, if there are zero modes at k = 0, they can be chosen 
to be eigenstates of V with eigenvalues ±1. At k = 0, 
for any eigenstate ipE,r of H(0) with energy E ^ 0, we 
can generate an eigenstate i^-e,t — V^}e,t with energy 
—E. Therefore, the only diagonal matrix elements of 
V in the basis of the eigenstates of Ti.(0) come from 
the zero-energy states. Let us denote by n± the num- 
ber of zero-energy states with P-eigenvalue ±1. Then, 
TtV = n + — n_. Therefore, TtV constitutes the lower 
bound on the number of zero energy states of Ho ■ 

Since TtV is independent of the basis, we can com- 
pute the trace in the coordinate basis. Consider again 
the square vortex lattice depicted in Fig. 1 with two 
hc/2e vortices per unit cell with the primitive vortex cell 
45° relative to the tight-binding lattice and the magnetic 
length I — 2(2n + 1). Then, there are only four spa- 
tial points which remain invariant under inversion and 
contribute to the trace. We listed the values of V in Ta- 
ble J] from which it follows that TtV = 4. Therefore, 
the index of 7i(0) is 4 and it has at least 4 zero-energy 
eigenstates. Same argument holds for HQs*) where 



k* G {(0,0),(7r/£,0),(7r/£,7r/^),(0,7r/£)}, with the only 
modification being Vv.* = (—l) x+y e lG r -f r Z^2, where 
G* = —2k*. So for I — 2(2n + 1) there are at least 
16 zero energy eigenstates in the first Brillouin zone. 

Consider now the same vortex arrangement, but for 
hc/e vortex lattice (or for hc/2e but ignoring the 
branch cuts). Then, j r = 1 everywhere and TrPk* = 
except at M point (k* = (n/t,Tr/£)) where Tr-p k . = 8, 
i.e. the index theorem guarantees only 8 zero energy 
states. This is the same number of zeros as in the ab- 
sence of magnetic field. Thus, while the effects of the 
superflow are perturbative, the effect of the branch cuts 
is non-perturbative. The symmetries discussed exist only 
if the bonds change sign exactly. Clearly, 16 zero energy 
states cannot be reached by perturbation theory around 
a state with 8 zero energy states, such as the <i-wave su- 
perconductor in the absence of the magnetic field. 

On the other hand, for £ = 2(2n), TrP = in both 
cases, because there are no lattice points invariant under 
inversion. The index theorem thus does not guarantee 
any zero energy states, i.e. the spectrum is gapped. 

Degenerate Dirac cones: Around each of the zeros in 
the magnetic Brillouin zone, the energy vanishes linearly 
with the deviation 5k from the degeneracy points. This 
follows from the standard k • p perturbation theory ar- 
guments as well as the fact that at /i = all the bands 
are doubly degenerate 0. Expanding the Hamiltonian 
7i(i5k) in the basis of the zero energy states near k = 
we find 

H ef f = V+Sk+Ti + VSk-Tt, (8) 

where r 1: 2 are 4x4 Dirac gamma matrices and dk± = 
(k x ± k v )/y/2. The spectrum near k* = thus consists 
of two doubly degenerate Dirac cones with energies 

£ = ±t^Vl5k 2 + + V 2 _5k 2 _. (9) 

The new nodal velocities V± must be determined from 
the wavefunctions of the degenerate multiplet. For the 
square vortex lattice in Fig^ we found analytically that 
for I = 2, Vi = 2(1 + a 2 ) ± -^(a 2 - 2(1 + ^2)), 
where the bare Dirac cone anisotropy a — Ao/t. Eq.© 
holds also near k* = (it/£,it/£), while near (w/1,0) and 
(0,n/£), Sk± — > ±5k T as required by the fourfold sym- 
metry. While £ — 2 corresponds to an unrealistically 
large magnetic field for physical parameters, £ = 2 is 
well representative of the sequence £ = An + 2 and has 
the virtue of being analytically soluble. We verified nu- 
merically that the qualitative results discussed below are 
indeed the same along the sequence £ = 4n+2, for integer 
n. It is not difficult to see that the effective Hamiltonian 
(JHJ) has the same form for any £ = An + 2. In addition, 
in the limit of large n, the Dirac velocities V± become 
universal functions of a in accordance with Simon-Lee 
scaling |4J. We verified this numerically (see FigJ^Jl. 
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FIG. 2: Upper panel: The energy dispersion in the gapless 
(nodal) (u = 0) and gapped (u 7^ 0) phases. Only the positive 
energies are shown as the spectrum is symmetric under E — > 
—E. For u = and £ — An + 2 each (red) band is doubly 
degenerate and there are 16 zero modes in the Brillouin zone. 
Finite fj, splits the degeneracy and the low energy spectrum 
corresponds to 8 massive Dirac fermions, each contributing 
equally to the spin Hall conductivity a% y — ±4^. 
Lower panel: The renormalized nodal velocities normalized to 
their value at I = 2 as functions of £ — An + 2, approach the 
universal value dependent on a = Ao/t only. Here a — 0.1. 

P-h symmetry breaking and topology of the gaps: 
For small but finite /j,, all the degeneracies are lifted 
Fig- ©1 but the spectrum is still symmetric under 
E — > —E for each k point due to the inversion symmetry 
of the vortex lattice Q . At low energies we have 

£ = ±\j£%+ml±\^^+ml (10) 

where the ± signs are uncorrelated, mo and mo are func- 
tions of £ and a, and £ 2 = t 2 (v\ 5k\ + Vl 5k 2 _ ). As 

shown in Fig.J5J), at each k* there is a local maximum in 
the lowest positive band of the spectrum; the local mini- 
mum is split and moved to {8k +1 Sk_} — {±fc mi „, 0} near 
(0,0) and (tt/1,tt/£), and {5k + ,5k_} = {0,±k min } near 

(n/£,0) and (0,n/£), where k mm = £JL^2|-±™£. 

In the vicinity of each of the 8 band-gap minima, the 
spectrum is equivalent to the 2+1 dimensional massive 
Dirac particle, with mass mo given by the gap minimum 
Emin- To see that, around each of the eight local energy 
minima, project onto the lowest energy basis, |±), such 
that H(k m i„)|±) = ±E m in\±), and expand the resulting 
effective Hamiltonian in 5k.: 

H e // = E min r 3 + v + 8k + Ti + vSk_T 2 , (11) 

Ti are the Pauli matrices acting on the 2-dimensional sub- 
space spanned by |±). Clearly, TL is equivalent to the 



Hamiltonian of a 2+1 dimensional massive Dirac particle. 
By the E — > —E symmetry of the spectrum, the nega- 
tive energy states are all occupied. It is well known, that 
in 2+1 dimensions, in the representation where v± > 0, 
the sign of the mass determines the Hall conductivity 
a xy — — isgn(mi))^- 0,0,0]- In the case at hand, the 
charge is not conserved due to broken U(l) symmetry of 
the superconductor, but the component of the spin along 
the magnetic field is conserved, and so spin Hall conduc- 
tivity, a* y , is quantized 0, H, • We verified that all 
8 of the Dirac particles have the same chirality, i.e. their 
additive contributions are the same. In order to deter- 
mine the total value of cr xy , generically, it is not enough 
to know the low energy part of the spectrum. Only the 
change of a s xy can thus be determined. Nevertheless, in 
this case, at fj, = 0, a xy — due to the particle- hole 
symmetry and since the change is ±8, a xy — ±4 (in 
appropriate units). We verified this numerically as well. 
The resulting state at \x ^ remains distinguishable from 
the case of hc/e vortex lattice, where \o~ xy \ < 2. 

While the above analysis pertains to the initial 
Ansatz for e l9rr ' , the main conclusions depend only 
on the translational and inversion invariance of A rr / x 
^2,rr'e~^^ r e _ 5^r' w hich we verified holds at each step of 
the self-consistent iteration for a model with the nearest 
neighbor attraction that stabilizes d-wave pairing. 

Finally, we conclude with comments about experi- 
ments. Recently it was observed that the low tem- 
perature thermal conductivity, k xx , of the ultraclean 
YBa2Cu307 showed pronounced deviations from the 
semiclassical predictions Rather, it was qualita- 

tively consistent with the onset of a new universal k, XXi 
distinct from its zero field value. As observed by Franz 
and one of us (O.V), the universal renormalization of the 
nodal velocities, shown in Fig[3 could provide a quali- 
tative explanation [20j . Presently, measurements of n xy 
extend to about 10K which is too high to observe the 
quantization discussed. 
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